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A novel approach to heteropolymer design is proposed. It is based on the criterion by Kurosky 
and Deutsch, with which the probability of a target conformation in a conformation space is 
maximized at low but finite temperature. The key feature of the proposed approach is the 
use of soft spins (fuzzy monomers) that leads to a design equation, which is an analog of the 
Boltzmann machine learning equation in the design problem. We implement an algorithm based 
on the design equation for the generalized HP model on the 3x3x3 cubic lattice and check its 
performance. 

KEYWORDS: design, heteropolymer, protein, optimization, learning, MTP criterion, Boltzmann machine, design 
equation, HP model 



§1. Introduction 

Recently, computer design of microscopic objects 
draws much attention of theoretical physicists. De- 
signing & helfimpplymer that fold into a given 
shape&BBlSMQ'&El' is the most challen ging one because 
it is a key technology for protein engineering. This prob- 
lem is also called "inverse folding" of heteropolymers. 
While the equilibrium conformation of a polymer with a 
given sequence of monomers is asked in a folding prob- 
lem, a sequence that folds into a given conformation of a 
polymer is requested in the corresponding inverse folding 
problem. 

The first step in the inverse folding problem is to for- 
mulate it-ias an optimization problem. In the pioneer- 
ing workEP of Shakhnovich and Gutin, a sequence <x of 
monomers (or "amino acids") that minimize the energy 
E(f\er) of a target conformation r is chosen under the 
condition of a fixed monomer composition. This crite- 
rion is based on the hypothesis that the energies of mis- 
folded conformations depend only on the composition of 
the polymer. Algorithms with this criterion have been 
shown to work in many practical problems with compu- 
tational efficiency. An optimal solution, however, is not 
always ensured ty 

Kurosky and Deutscbdllf proposed a different crite- 
rion. In their framework,E2P the equilibrium probability 

Z p =J2eM-PE(r\<r)) (1.2) 
r 
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that a polymer is found in a target conformation f is 
maximized at a sufficiently low temperature 1//3. Here 
y^ r means the summation over all possible conforma- 
tions. Hereafter, we denote this criterion as "MTP cri- 
terion" (Maximum Target Probability criterion). The 
MTP criterion automatically excludes solutions with de- 
generate ground states when there exists other solution 
with an unique ground state. In addition, if we keep the 
temperature 1//3 low but a finite value, it selects a so- 
lution with a larger energy gap above the ground state. 
An advantage of MTP criterion is that it is useful to sys- 
tematically understand various algorithms as methods 
for approximate maximization of MTP criterion, for ex- 
ample, the original algorithm by Deutsch and r Kurosky,o' 
the algorithm by Morrissey and ShakhnovichcP , and the 
algorithms by Seno et alu 1 

The purpose of this paper is to present a new approach 
to the inverse folding problem based on the MTP crite- 
rion. At the present status, we are mostly concentrated 
on the design of the sequences that has a given target 
conformation as the unique ground state. For this pur- 
pose, the proposed approach has a considerable advan- 
tage over other approaches, e.g., the one in ref. ^, which 
will be described in the section || and ^. Although we 
are interested in the low temperature limit /? — > oo, the 
use of MTP criterion at low but finite temperature is 
essential in the derivation of the proposed algorithm. 

A goal of inverse folding of hetropolymers is to con- 
trol the dynamical property of the designed sequence as 
well as thermodynamic one. That is, we want to design 
a sequence that fold fast into the desired conformation. 
Around this problem, there are a number of issues on 
the relation between the thermodynamic property and 
dynamical property of hctcropolymecs .-below the fold- 
ing temperature: While some authorstilEP argue that the 
foldability of a heteropolymer mostly depends on the en- 
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ergy gap between the ground state and lower excited 
states, the other authorala' claim the relevance of some 
different properties. In our context, it will be interesting 
to test the dynamical properties of the sequence designed 
by the proposed algorithm. Wc will, however, leave it for 
a future problem and restrict ourselves to an innovation 
of the computational technique in this paper. We also 
hope that the development of an efficient algorithm will 
contribute investigations into the above-mentioned prob- 
lem. 

§2. The Design Equation Approach 

Let us discuss the proposed approach - the design 
equation approach. Our starting point is the analogy 
to the learning in artificial neural netw orks . In fact, 
the maximization of the probabilitY-ecL ( |l . l| ) is analo- 
gous to Boltzmann machine learningcJ'O in neural net- 
work theories (and also maximum likelihood estimation 
in statisticsEj) ) , where the parameters that maximize the 
probability of a given data are selected as an optimal so- 
lution. Although this analogy to Boltzmann machine 
learning was already mentioned in refs. ^ and ||, it has 
not been fully explored in the literature. Specifically, no 
direct analog of Boltzmann machine learning equation 
has been discussed in connection with the inverse folding 
problem. Here we further pursue the analogy between 
the MTP criterion and Boltzmann machine learning and 
give a new approach based on a design equation, an ana- 
log of Boltzmann machine learning equation. For this 
purpose, we introduce "soft spins" or "fuzzy monomers" 
as a tool of computation. 

We will discuss pur .approach in terms of the general- 
ized HP mode£3'L3'E3EH' of protein. This model consists 
of a self-avoiding polymer chain on a lattice with two 
types of monomers indexed by 1 and —I. As a model 
of protein, the indices ±1 indicate hydrophobic (H) and 
polar (P) residues, respectively. The interaction energy 
U between a pair of monomers (the contact energy) is 
defined as 17(1,1) = ei, 17(1, — 1) = Z7(— 1, 1) = e 2 , and 
U(— 1, — 1) = C3. It acts only between monomers on near- 
est neighbor sites but not consecutive along the chain. 
The energy E(r\er) of a conformation r = {r{\ (r"j rep- 
resents the vector of coordinates of the ith monomer) of 
a polymer with a sequence <t — {o~i} is written as 



E{r\cT) = \Y. U ^ (T i) 6 ^- r >\- 1 ) 



Vij 



(2.1) 



where er^ € { — 1,1} represents the type of ith monomer. 
The factor rjij G {0, 1} takes the value if and only if 
monomers i and j are consecutive along the chain. For 
later convenience, we note that the interaction energy 
U(o~i, o~j) of the model is written as 



U((Ti,0-j 



ei - 2e 2 + £3 



et >( 



(2.2) 



ei - £3 ei + 2e 2 + £3 
— ^ — > £ c = : 



(2.3) 

At this point, we introduce "soft spins" m = {to^} 
each of that takes a continuous value — 1 < to, < 1 and 



substitute the original binary variables <x = {(7^}. Non- 
integer values of the variable m, ( "fuzzy monomer" ) have 
no physical meaning in the generalized HP model. They 
are, however, convenient tools for the-jComputation, as 
they are in the Hopfield-Tank methodEP for combinato- 
rial optimization problems. We take the following form 
eq. (^4|) of the energy function, which is a straightfor- 
ward extension of eq. (2.1) to the soft spin model. 



E*(r\m) = ^^2u m (m i ,m j )6{\r i 



1)% (2-4) 



The only constraint imposed on the modified energy 
function E*(r\m) with continuous variables m = {to^} 
is that it coincides with the original energy function 
cq. ( |2.l| ) when |mj| = 1 for all i. In this paper, we 
consider two possible choices for U* , 

e a mimj + ^(m 4 + rrij) + e c t (2.5) 



and 

U 2 (mi,mj) = 



e a rrnmj 



I m j I ) 



— (riii |to.j 4 ,1 , ! in, 

(2.6) 

is an obvious extension 



The former expression eq. (|2.5[ 
of the original energy function eq. ([2.2| ) and it is easy to 
see that the latter expression eq. ( |2.6| ) is also coincides 
with the original energy function when \m,i\ = I to, I = 1. 



If we substitute these energy functions in eq. (1.1) of 
the MTP criterion, it gives 



Pp(r\m) 



Z B = 



r 



exp(-PE*(r\m)) 



exp(-/iE*(r|m)) 



(2.7) 



(2.8) 



The variables {to^} that maximize the expression 
cq. ( |2.7| ) take, in general, non-integer values and not 
necessary satisfy the relations — 1 < < 1. Then we 
introduce a penalty term 

A x ^ 



V p {m) 



4^ 



{mi - If 



(2.9) 



to force them towards 1 or — 1. The value of the param- 
eter A controls the strength of the penalty. Using this 
penalty term, we arrive at the cost function 



V cost {m) = - log Pp(r\m) + V p (m) 



(2.10) 



to be minimized in the present approach. 

The use of the soft spins and the penalty term are 
key features of our approach. By virtue of them, we can 
differentiate the cost function (2.10) with to,; and write 
down a set of equations, 

drrii dVcost 



dt 



dm; 



U<fi,m) - \mi(mi - 1) , (2.11) 



/ i (/3,m) = /3^ 
j 

x {5{\ri -r-j 



dll*(mi, Tfij ) 
dmi 

1) - (S(\ri 



■ rj 1-1))/,} (2.12) 
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that minimize the cost function eq. (2.10) with a gradient 
decent method. Here the variable t is a fictitious time 
and the constant r controls the time scale. The average 



-l))/3 indicates the canonical average of S(\ri 



1) at the inverse temperature /3, i.e. 



i)> 







£ 



S(\n- 



l)Pp(r\m) . (2.13) 



We omit the factors rjij in the expression of fi in 
eq. (2.12) because their effects cancel between the first 
and second terms in the brace { } in eq. (2.12). 

In this paper, we refer to the set of equations, 
eqs. ( 2.11 ) and (2.12), as the design equation for this 



problem. The design equation is nothing but an analog 
in the inverse folding problem of the Boltzmann machine 
learning equation. When the value of the control param- 
eter A is gradually increased to +oo as the fictitious time 
t — > oo, the value of each soft spin rrii converges to ±1, 
which defines a sequence with proper meaning in the 
original problem. It is easy to see that a sequence whose 
unique ground state coincides with the desired confor- 
mation satisfies the equation dV cos t/dmi = 0. Thus, 
for sufficiently low temperature 1/(3, the output of the 
present procedure is a candidate for the solution of the 
original problem with the MTP criterion. 

Details of our implementation of the design equation 
are shown in the followings: 

1. Initialization. 



Set k := 1, 



A := A . 



2. Calculation of the canonical averages. 
Calculate {8{\ri — r,-| — by the exact enumera- 
tion or by a Monte Carlo procedure. The former is 
possible only for short chains. 

3. An iteration of the discretized design equation. 



dU*(rrii 1 mj) 



x{«(|r i -r i |-l)-(% i -rj|-l)> /s } 
m l := m l + A {/,, - Xm^mf - 1)} . 

4. Clipping. 

If mi > 1, set rrii := 1. 
If mi < —1, set rrii '■= — 1- 

5. Increment the variables and check the convergence. 
Set k := k + 1 and A := A + a. If a prescribed stop- 
ping criterion is satisfied or the counter k exceeds 
a prescribed maximum of the iteration, end up the 
calculation. Else, return to step 2. 

The constant a controls the increase of the strength A of 
the penalty and the constant A controls the size of steps 
in the gradient decent. The tunable control parameters 
are Ao, a, A, {3 and initial conditions {m®}. 

There are two candidates of the stopping criterion. A 
possible criterion is "If, for all i, the value of m, is suf- 
ficiently close to 1 or —1 and does not change in several 
consecutive iterations, end up the calculation" . If we can 
check at every step whether the output is a solution of 
the problem, we can use another criterion based on the 



"forced discretization" ai of defined as follows: 

(jj =1 if nii > 0, 

<7j = — 1 otherwise. (2-14) 

With this stopping criterion, we end up the calculation 
before reaching prescribed maximum number of the iter- 
ation if and only if the target conformation is the unique 
ground state of the polymer with the sequence . 

§3. Related Works 

In this section, we discuss some of earlier works in con- 
nection with the design equation approach. First, the 
use of the soft spin variables is, in itself, not a new tech- 
nique for the protein design. They are also introduced 
as optimization tools in the refs. ||,^l],||. In these refer- 
ences, however, criteria similar to that of Shakhnovich 
and Gutirff are used and the MTP criterion is not em- 
ployed. The design equation is not appeared in these 
works. 

On the other hand, there have been proposed sev- 
eral algorithms based on the MTP criterion, which use 
different approximations and computational techniques. 
In the original algorithm proposed by JDeutsch and 
KuroskyeP (see also Kurosky and DeutscbEP), the loga- 
rithm of the partition function Zp is approximated by the 
first order cumulant expansion and equally-weighted av- 
erage over all compact self-avoiding conformations of the 
polymer. The resultant expression is optimized through 
a simulated annealing procedure. For the generalized 
HP model, the high temperature approximation leads to 
the expression apparently similar to the the right hand 
side of the design equation. There is, however, a major 
difference, because the right hand side of the design equa- 
tion eq. (2.12) can automatically incorporate the effect 
of higher order correlations in the conformation space as 
the magnitudes of {|?rii|} are increased. Morrissey and 
ShakhnovicbE' also gives an algorithm based on MTP 
criterion, which uses a higher order cumulant expansion 
of the free energy. These cumulants are evaluated by 
a mean field approximation in conformation space, so 
that it still keeps the computational economy of the al- 
gorithms based on the Shakhnovich and Gutin's crite- 
rion™ Its performance seems mostly dependent on the 
validity of the approximation in the conformation space. 
We will discuss another aspect of their work at the end 
of the paper.- 

Seno et al& developed a dual Monte Carlo algorithm 
that is most faithful to the MTP criterion. In their al- 
gorithm the partition function Zp is calculated by an 
important sampling in the conformation space and the 
calculated value of P(r\cr) is optimized through a simu- 
lated annealing procedure in the sequence space. They 
also test the performance of the algorithm where the 
Monte Carlo calculation of the partition function Zp is 
replaced by the exact enumeration of conformations. In 
the next section, we compare the performance of an al- 
gorithm based on the design equation with this version 
of algorithm. A significant difference between our ap- 
proach and theirs is that the ave rag e (S(\ri — rj\ — l))p 
over the Gibbs distribution eq. ( |2.7| ) is required in our 



Author Name 



approach, instead of the partition function Zp required 
in the algorithms by Seno et al. 

For any algorithm based on the simulated annealing in 
the sequence space, the value of a cost function should be 
evaluated once per a spin flip in the simulated annealing 
procedure. For example, in the algorithms of Seno et 
al, the value of the partition function Zp is calculated 
once per a trial of changing the type of one monomer 
in the sequence. In the design equation approach, the 
canonical averages (S(\ri — rj\ — 1))^ is calculated only 
once per an iteration of the discretized design equation. 
As will be shown later, this difference causes a significant 
advantage of the algorithms based on the design equation 
when the computational cost of the cost function and its 
derivatives are comparable and both intensive. 

§4. Numerical Experiments 

Here, we test the design equation approach for a three 
dimensional generalized HP model ( 2T ) on the cubic lat- 
tice with 27 monomers. In the experiments, we restrict 
the conformations of the polymer to maximally compact 
self-avoiding ones filling the 3x3x3 lattice. Numerical 
experiments with compact conformations on small lat- 
tices are commnnJn the study of protein foldingE2P and 
inverse foldingptfla) because then the exact enumeration 
of the conformations is possible. 

In the following experiments, the canonical averages 
(S(\ri — Tj \ — l))/3 are calculated by the exact enumera- 
tion. The parameters (£1,62) = (— 1,0) and (—2.3,-1) 
are used. We set £3 = throughout the experiments. 
Note that if we restrict .our attention to the maximally 
compact conformations,!^ the result is not affected by 
the addition of a common constant to all of the parame- 
ters £1, £2, £3 (equivalently, by the value of the parameter 
£ c ). Thus, we can set £3 = without any loss of gener- 
ality. We test both of the modi fied energy function U 1 , 
U 2 defined in eq. (2.t) and eq. (2.5). 

In Fig.l, an example of the designing process is shown. 
In this case, the energy function U 1 is used and the algo- 
rithm successfully found a sequence that has the target 
conformation as the unique ground state (a "good" se- 
quence) . 

To check the performance of the present algorithm, 
we execute the test calculations for five different target 
conformations in Table I. The conformation #1 is "the 
most designable one" for the generalized HP model on 
the 3x3x3 lattice according to the definition in ref.|||. It 
seems to be one of the easiest target also in the present 
cases. The conformations #2 - #5 are randomly cho- 
sen ones and some of them are much more difficult as a 
target than #1. For each target conformation, 200 trials 
with different initial sequences {m? } are performed. The 
number of the successful attempts that result in a "good" 
sequence is recorded. We also record the number of the 
iterations needed to reach a good sequence when we use 
"forced dicretization" as a stopping criterion. The ini- 
tial values {m®} of {mi} are generated randomly within 
the range [—10,10] using a uniform random number. The 
temperature 1/(3 in MTP criterion is set to 0.01. The 
values of the control parameters in the design equation 
are A = 0.5, A = 0.0, a = 0.5. 



To get a feeling of the difficulty of the problem, we 
also implement an algorithm based on simulated anneal- 
ing in the sequence space. This algorithm, to which we 
will refer as "SA" , is essentially an algorithm by Seno et 
al.U with the calculation of Zp by the exact enumera- 
tion. In the experiments of SA, 10 trials with different 
initial sequences are performed for each target confor- 
mation. The fictitious temperature Tsa of the simulated 
annealing is initially set as Tg A = 2.0 and updated by 
Tsa ■— 0.8 x Tsa in every Monte Carlo step per spin. 
Note that the fictitious temperature Tsa, which is in- 
troduced for the optimization in the sequence space, is 
independent of the tem pera ture 1/(3 that determines the 
Gibbs distribution eq. (1.1) in the conformation space. 
The value of 1/(3 is kept 0.01 throughout the experi- 
ments. 

The results of the experiments with (£i,£2) = (—1,0) 
and (—2.3, —1) are shown in Table II and Table III, re- 
spectably. 



#1 


brf2 1 2U2 i>2 T2 h ldrb 2 h dfuru 


#2 


r 2 f2 ubufl2 d,2 brfubulbr2 dl%j 


#3 


r 2 f2ublulfd 2 br fu2rb2dluldf 


#4 


r 2 f2U-l2ub2druf 2 rb2dfldflbu 


#5 


r 2 f2ubldfu2rblbrdl2uf2d2 bu 



Table I. The target conformations in the experiments. Each con- 
formation is encoded by the symbols r(ight), i(eft), u(p), d(own), 
/(forward), fc(ackward). The expressions ri, f2 ■■ are contracted 
forms of rr,ff etc. 



The design equation with the modified energy func- 
tion U 2 successfully finds at least one good sequence for 
each of the five targets both for (fi,£2) = (—1,0) and 
(—2.3,-1). The comparison to SA shows that the cost 
S of the calculation in each successful run is much less in 
the proposed algorithm than in SA. This result is highly 
dependent on the computational advantage of the design 
equation approach, i.e, the economy of the evaluation of 
the cost function that we have already mentioned at the 
end of the previous section. Note that this advantage 
will be larger as the length of the polymer increases. On 
the other hand, the rate R of the success is lower in the 
proposed algorithm in most cases. Then overall perfor- 
mance seems comparable in both algorithms. If we define 
an index P = 100 • R/S as a measure of efficiency, the 
proposed algorithm with U 2 makes better scores than 
SA in most cases. 

Our experiments show that the performance with 
the energy function U 1 is rather poor in the case of 
( e i,£2) = ( — 2.3,—l). On the other hand, the algorithm 
works well with the energy function U 2 in both cases 
(ei, e 2 ) = (—1,0) and (—2.3,-1). A possible reason of 
this difference is that the energy U 2 is a better repre- 
sentation of the original energy function eq. (2.2) in the 
sense that the factor miirij in the first term and the fac- 
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Fig. 1. An example of the designing process. 

The value of each mi is plotted versus the number of the iteration. A successful convergence to a good sequence of the target 
conformation #1 (see Table I) is shown. Note that "forced discretization" provides a good sequence with 20 iterations of the discretized 
design equation. Here (ei,£2) = ( — I, 0) and A = 0.5, Aq = 0.0, a = 0.5. 



tor mi|mj| + mj\rrij\ in the second term of eq. (2.6) are 
both in the same order of magnitude for small values 
of |toj|s. In the en ergy function U , the second term 
y(mi + nij) in eq. (2.5) dominates when the values of 
|m,|s are small and the early stage of the designing pro- 
cess is mostly driven by the second term when we start 
from an initial condition with small |mj|s. 

§5. Summary and Future Problems 

In Summary, we have proposed a novel approach for 
hetcropolymer design that is based on the Maximum 
Target Probability criterion by Kurosky and Deutsch. 
The essential point of the proposed approach is the in- 
troduction of soft spins (fuzzy monomers) that leads to 
a design equation, which is an analog of the Boltzmann 
machine learning equation in heteropolymer design. We 
have tested this approach, which we call the design equa- 
tion approach, for a generalized HP model on the 3x3x3 
lattice and have shown that it could successfully find 
good sequences for target conformations with different 
degrees of the difficulty. With these examples, our im- 
plementation of the design equation approach shows at 



least comparable performance to an algorithm based on 
simulated annealing in sequence space, when a suitable 
form of the modified energy function U* is chosen. 

In this paper, we have evaluated the performace of 
the algorithm by the ability to give a "good sequence" , 
a sequence whose unique ground state coincides with a 
desired connformation. As has been mentioned in the 
introduction of the paper, it is also interesting to test dy- 
namical property of the outputs of the algorithm below 
their-folding temperatures. On the other hand, a recent 
workEr of Morrissey and Shakhnovich arises a problem 
of selecting sequences whose equilibrium state at a fi- 
nite, not necessarily low, temperature is dominated by a 
given shape. The design equation approach might also 
be useful in this problem. For this purpose, however, 
further study of the behavior of the design equation at 
finite temperature 1//3 is required. 

The design equation approach is fairly general and var- 
ious modifications are possible. Here we touch on a few 
important extensions. First, we can substitute a Monte 
Carlo simulation for the exact enumeration in the present 
algorithm. This is important because exact enumeration 
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Target 




U 1 






u 2 






SA 






R 


S 


P 


R 


s 


P 


R 


S 


P 


#1 


0.82W 


12.2 


6.7 


0.56(t) 


16.4 


3.4 


1.0 


162.0 


0.62 


#2 


0.215M 


15.0 


1.4 


0.16<*) 


14.8 


1.1 


0.9 


294.3 


0.31 




0.20<t) 


18.5 


1.1 


0.25(t) 


14.0 


1.8 








#3 


0.30<*) 


16.0 


1.9 


0.355(t) 


15.9 


2.2 


0.8 


162.0 


0.49 


#4 


0.04<*) 


17.6 


0.23 


0.055(t) 


14.4 


0.38 


0.4 


256.5 


0.16 


#5 


0.08M 


12.1 


0.66 


o.uct) 


14.7 


0.75 


0.1 


232.2 


0.043 



Table II. The results with (c i , oj) = (—1,0) for the target conformations in Table I. Rates (R) of finding good sequences, average 
numbers (5) of iterations (or spin flips in SA) and "efficiency" P = 100 ■ R/S are shown in the cases of (1) the proposed algorithm with 
the energy function U 1 , (2) the proposed algorithm with the energy function U 2 , and (3) an algorithm based on simulated annealing 
in the sequence space (SA). In SA, total number of spin flips until a "good" sequence is used as a correspondence of the number of 
the iteration of the design equation, because the calculation of the partition function is required in each spin flip in SA. The initial 
condition is generated randomly within (— w,w), where (*) w = 1.0, (f) w = 0.1. 



Target 




U 1 






U 2 






SA 






R 


s 


P 


R 


S 


P 


R 


S 


P 


#1 


LOW 


8.5 


11.8 


0.75(t) 


10.0 


7.5 


1.0 


105.3 


0.95 


#2 


0.005M 


1.0 


0.5 


0.005M 


12.0 


0.04 


0.9 


186.3 


0.48 




o.oi(t) 


2.5 


0.4 


o.i(t) 


10.8 


0.93 








#3 


0.002M 


5.0 


0.04 


0.225<t) 


13.0 


1.7 


0.8 


210.6 


0.38 


#4 


O.OOM 






0.015<t) 


4.8 


0.31 


0.0 






#5 


0.00<*> 






0.045<t) 


11.7 


0.38 


0.3 


251.3 


0.12 



Table III. The results with (ei,e2) = (—2.3, —1) . The set of the targets is the same as that in the Table II (see Table I). The meaning 
of the notations are shown in the caption of Table II. 



of the conformations is impossible for longer or off-lattice 
polymers. Research in this direction is now in progress 
and we can successfully implement an algorithm based on 
the design equation with a dynamical Monte Carlo sim- 
ulation—The results will be reported in the forthcoming 
paper .Ef 

Another important challenge is the extension to the 
cases with monomers of many different types (many let- 
ter cases). Although a formal extension to many letter 
cases is not difficult, the test of the performance of the 
design equation in such cases is also left for the future 
study. 
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